Multi-arm Spacecraft Model Predictive Control Method Based on Mixture of Gaussian Processes, Equipment, and Medium

ABSTRACT

The present disclosure provides a multi-arm spacecraft model predictive control method based on the mixture of Gaussian processes, equipment, and a medium. Model predictive control has excellent performance in dealing with complex nonlinear systems such as multi-arm spacecrafts with various constraints, and is widely applied to ground robots, unmanned aerial vehicles, autonomous driving and other practical scenarios. Therefore, a task space controller is designed based on the model predictive control in the present disclosure. Besides, in order to enhance the anti-interference capability of the present disclosure, an interference model is established and compensation is carried out in the model predictive control by utilizing the characteristics of small training data volume and high training speed in the mixture of Gaussian processes. Finally, a thrust distribution method is designed to complete platform control. The method provided by the present disclosure is convenient and intuitive in design and has relatively high practicability.

TECHNICAL FIELD

The present disclosure belongs to the technical fields of on-orbitservices, spacecraft control, and model predictive control, and inparticular relates to a multi-arm spacecraft model predictive controlmethod based on the mixture of Gaussian processes, equipment, and amedium.

BACKGROUND OF THE PRESENT DISCLOSURE

In recent years, an on-orbit service technology has received widespreadattention from countries around the world. Major space agencies such asChina, the United States, ESA, and Japan have conducted a large numberof ground and space experiments, achieving various types of on-orbitservice tasks such as satellite capture, refueling, componentreplacement, and space debris cleanup. However, the spacecraft on-orbitservice tasks are all described in a task/Cartesian space, while thetraditional controller design is mostly expressed in a joint space,which requires additional inverse kinematics solution, so that besidesthe lack of intuitiveness, there are problems such as mismatch ofcontrol ability and large amount of calculation. Therefore, the presentdisclosure proposes a task space controller design method based on modelpredictive control (MPC), which is simple, direct, and suitable forprocessing high-dimensional nonlinear models of a multi-arm spacecraft.

Furthermore, there are various disturbances in the process of conductingexperiments on the ground and executing tasks in space for an on-orbitservice spacecraft, including parameter uncertainty, spatial andenvironmental torque disturbances, and input and measurement noise,which have a significant impact on the performance of on-orbit servicetasks. Traditional solutions have their own advantages anddisadvantages: the method based on a robust design needs to assume thatthe total disturbance has an upper bound; the method based on anadaptive disturbance observer is capable of estimating disturbancesonline, but it has higher requirements for initial values and is poor inperformance at initial moments; and the method for constructingdisturbance models based on neural networks and fuzzy networks requiresoffline training and is better in initial performance, but it requires alarge amount of training data. Considering the actual workflow of theon-orbit service spacecraft, it is promising to design an off-linedisturbance model construction method with small data demand and fasttraining speed by using the Mixture of Gaussian Processes (MGP).

Finally, the traditional on-orbit service spacecrafts are generally in afree-floating mode, which do not control the satellite platform, or onlyrely on flywheels for platform attitude adjustment, resulting in weakmaneuverability and limited operating range, which affects theefficiency and performance of on-orbit services. Aiming at this problem,the present disclosure employs jet thrusters to adjust the spacecraftplatform pose and designs a corresponding thruster control distributionmethod.

SUMMARY OF THE PRESENT DISCLOSURE

The purpose of the present disclosure is to propose a multi-armspacecraft model predictive control method based on the mixture ofGaussian processes, equipment and a medium for the deficiencies of theexisting multi-arm spacecraft control methods.

The present disclosure is realized by the following technical solution,and the present disclosure proposes a multi-arm spacecraft modelpredictive control method based on the mixture of Gaussian processes,which specifically includes:

-   -   establishing dynamic and kinematic models of a multi-arm        spacecraft as a prediction model in MPC;    -   for an estimated value {circumflex over (d)}of a disturbance        term in the prediction model, inputting initial excitation to        establish a data set, and employing the mixture of Gaussian        processes for training to obtain residual dynamics of the        multi-arm spacecraft as the estimated value of the disturbance        term;    -   according to actual drive input saturation constraints,        constructing a model predictive controller to implement the        tracking of a desired trajectory for end poses and platform pose        of multi-arm spacecraft manipulators, where the saturation        includes saturation of joint motors and thrust saturation of        platform thrusters are considered; and    -   setting the platform thrusters to be in an on-off control mode,        and assigning a platform continuous control instruction        generated by the MPC as start-up time of each of the thrusters        to obtain a final thruster driving instruction.

Further, for a multi-arm spacecraft with n manipulators, each of whichhas m degrees of freedom, and with a platform equipped with l jetthrusters, the dynamic model of the multi-arm spacecraft is firstdescribed in a joint space:

-   -   the dynamic model of the multi-arm spacecraft is expressed as:

M{umlaut over (Q)}+C+d=u   (1)

where M represents an inertial parameter of the multi-arm spacecraft; Crepresents a nonlinear term; d is the disturbance term; Q=[r₀ ^(T),q^(T)]^(T) represents information about the platform pose and jointangle, r₀ is the 6-degree-of-freedom platform pose, and q=[q₁, . . . ,q_(n)]^(T) is the joint angle; and u=[F₀ ^(T), F_(q) ^(T)]^(T)represents a control force and a control torque, F₀ represents aplatform control force/torque, and F_(q) represents a joint torque.

Further, based on a mapping relationship between the end-effectors ofthe manipulators and the states of joints and a platform, thedescription of the end-effector poses of the on-orbit service spacecraftis obtained by using a DH parameter method:

$\begin{matrix}{\begin{bmatrix}r_{0} \\r_{e1} \\ \vdots \\r_{en}\end{bmatrix} = {\begin{bmatrix}r_{0} \\{{\,_{}^{1}T_{0}^{1}}{\,_{}^{1}T_{1}^{2}}\ldots{\,_{}^{1}T_{m - 1}^{m}}} \\ \vdots \\{{\,_{}^{n}T_{0}^{1}}{\,_{}^{n}T_{1}^{2}}\ldots{\,_{}^{n}T_{m - 1}^{m}}}\end{bmatrix} = {g(Q)}}} & (2)\end{matrix}$

where a homogeneous transformation matrix ^(n)T_(m-1) ^(m) is definedbased on DH parameters θ_(m),d_(m),a_(m),α_(m):

$\begin{matrix}{{\,_{}^{n}T_{m - 1}^{m}} = \begin{bmatrix}{\cos\theta_{m}} & {{- \sin}\theta_{m}} & 0 & a_{m} \\{\sin\theta_{m}\cos\alpha_{m}} & {\cos\theta_{m}\cos\alpha_{m}} & {{- \sin}\alpha_{m}} & {{- d_{i}}\sin\alpha_{m}} \\{\sin\theta_{m}\sin\alpha_{m}} & {\cos\theta_{m}\sin\alpha_{m}} & {\cos\alpha_{m}} & {d_{i}\cos\alpha_{m}} \\0 & 0 & 0 & 1\end{bmatrix}} & (3)\end{matrix}$

finally, a state space expression is built:

$\begin{matrix}\left\{ \begin{matrix}{\overset{˙}{X} = {{f\left( {X,u} \right)} = \begin{bmatrix}\overset{.}{Q} \\{M^{- 1}\left( {\tau - C - \hat{d}} \right)}\end{bmatrix}}} \\{\xi = {g(Q)}}\end{matrix} \right. & (4)\end{matrix}$

where X=[Q^(T), {dot over (Q)}^(T)]^(T), ξ=[r₀ ^(T), r_(el) ^(T), . . ., r_(en) ^(T)]^(T) represents the platform pose and the end-effectorposes of the manipulators, and {circumflex over (d)} represents theestimated value of the disturbance term; and

-   -   sampling time Δt is set, the continuous state space expression        is converted into a discrete time model, and system states at a        future moment k+N are predicted based on a value measured at a        current moment k:

$\begin{matrix}\left\{ {\begin{matrix}{X_{k + 1} = {X_{k} + {\Delta{t \cdot {f\left( {X_{k},u_{k}} \right)}}}}} \\{\xi_{k} = {g\left( Q_{k} \right)}}\end{matrix}.} \right. & (5)\end{matrix}$

Further, the recorded p groups of data sets are used for training of themixture of Gaussian processes, where an input x contains the multi-armspacecraft pose Q, speed {dot over (Q)}, and the increment Δu of acontrol variable, and an output y is a difference value between anactual angular velocity and an angular acceleration calculated by anominal model, that is, the residual dynamics, which is used to capturethe impact {circumflex over (d)} of the total disturbance on the system;

-   -   the Gaussian processes assume in advance that samples obey a        Gaussian distribution:

y(x)˜N (μ(x), k (x,x′))   (6)

where μ (x) represents a mean function; k (x, x′) represents acovariance function, and an expression thereof is as follows:

$\begin{matrix}{{k\left( {x,x^{\prime}} \right)} = {{\sigma_{f}^{2}{\exp\left\lbrack \frac{- \left( {x - x^{\prime}} \right)^{2}}{2l^{2}} \right\rbrack}} + {\sigma_{n}^{2}{\delta\left( {x - x^{\prime}} \right)}}}} & (7)\end{matrix}$

where σ_(n) is equal to 0.3, σ_(f) and l are hyperparameters to besolved, and

${\delta\left( {x - x^{\prime}} \right)} = \left\{ {\begin{matrix}{1,} & {x = x^{\prime}} \\{0,} & {x \neq x^{\prime}}\end{matrix};} \right.$

the covariance function is solved by using an optimal maximum likelihoodestimation method:

$\begin{matrix}{{\max\limits_{\sigma_{f},l}\log\left( {y{❘x}} \right)} = {{{- \frac{1}{2}}y^{T}K^{- 1}y} - {\frac{1}{2}\log{❘K❘}} - {\frac{n}{2}\log\left( {2\pi} \right)}}} & (8)\end{matrix}$

after the hyperparameters of the Gaussian processes are obtained viatraining, the corresponding output value y_(*) can be solved by usingthe input information x_(*) of a prediction point; a joint probabilitydensity function of the samples and a test set is firstly determined:

$\begin{matrix}{\begin{bmatrix}{y(x)} \\y_{*}\end{bmatrix} \sim {N\left( {0,\begin{bmatrix}K & K_{*}^{T} \\K_{\star} & K_{\star *}\end{bmatrix}} \right)}} & (9)\end{matrix}$ where $\begin{matrix}{K_{**} = {k\left( {x_{*},x_{*}} \right)}} & (10)\end{matrix}$ $\begin{matrix}{K_{*} = \left\lbrack \begin{matrix}{k\left( {x_{*},x_{1}} \right)} & {{k\left( {x_{*},x_{2}} \right)}\ } & \ldots & \left. {k\left( {x_{*},x_{p}} \right)} \right\rbrack\end{matrix}\  \right.} & (11)\end{matrix}$ $\begin{matrix}{K = \begin{bmatrix}{k\left( {x_{1},x_{1}} \right)} & {k\left( {x_{1},x_{2}} \right)} & \ldots & {k\left( {x_{1},x_{p}} \right)} \\{k\left( {x_{2},x_{1}} \right)} & {k\left( {x_{2},x_{2}} \right)} & \ldots & {k\left( {x_{2},x_{p}} \right)} \\ \vdots & \vdots & \ddots & \vdots \\{k\left( {x_{p},x_{1}} \right)} & {k\left( {x_{p},x_{2}} \right)} & \ldots & {k\left( {x_{p},x_{p}} \right)}\end{bmatrix}} & (12)\end{matrix}$

-   -   a posterior probability expression of the prediction point is        obtained based on Bayesian reasoning:

y _(*) |y(x)˜N(K _(*) K ⁻¹ y, K _(**) −K _(*) K ⁻¹ K _(*) ^(T))   (13)

a predicted value of the Gaussian processes is accordingly obtained:

y _(*) =K _(*) K ⁻¹ y   (14).

Further, a local Gaussian process method is utilized in the Gaussianprocesses of the multi-arm spacecraft. The local Gaussian process methodincludes: clustering based on a Gaussian mixture model method,performing regression on each of Gaussian components, and finally fusinglocal predicted values during prediction;

-   -   assuming that the input x obeys the mixture distribution of the        M Gaussian components, a probability density distribution        function of the input x is as follows:

$\begin{matrix}{{p\left( {x{❘\Theta}} \right)} = {\sum\limits_{i = 1}^{M}{\pi_{i}{p\left( {x{❘{\mu_{i},\sum_{i}}}} \right)}}}} & (15)\end{matrix}$

where π_(i) is a mixing coefficient, which represents the probabilitythat the sample comes from the i^(th) Gaussian component; p(x|μ_(i),Σ_(i)) represents the probability that the sample is generated by thei^(th) Gaussian component; and {μ_(i), Σ_(i)} is the hyperparameter ofthe mixture of Gaussian processes.

Further, after the training samples are obtained, the model parametersare first initialized, and then the model parameters are updated bymeans of the continuous iteration of steps E and M ofexpectation-maximization until the model converges,

-   -   where step E uses the parameters {π_(i), μ_(i), Σ_(i)} estimated        in step M to calculate a posterior probability r_(i) of the        sample coming from the i^(th) Gaussian component:

r i = π i ⁢ p ⁡ ( x ⁢ ❘ "\[LeftBracketingBar]" μ i , ∑ i ) ∑ j = 1 M π j ⁢ p⁡( x ⁢ ❘ "\[LeftBracketingBar]" μ j , ∑ j ) ( 16 )

in step M, the posterior probability r_(i) of the training samples iscalculated through step E, and the model parameters are updated bymaximizing a logarithmic likelihood function:

$\begin{matrix}\left\{ \begin{matrix}{\mu_{i} = {\sum\limits_{k = 1}^{n}{r_{k}x_{k}/{\sum\limits_{k = 1}^{n}r_{k}}}}} \\{\Sigma_{i} = {\sum\limits_{k = 1}^{n}{{r_{k}\left( {x_{k} - \mu_{i}} \right)}\left( {x_{k} - \mu_{i}} \right)^{T}/{\sum\limits_{k = 1}^{n}r_{k}}}}} \\{\pi_{i} = {\frac{1}{n}{\sum\limits_{k = 1}^{n}r_{k}}}}\end{matrix} \right. & (17)\end{matrix}$

for each group of Gaussian components obtained by clustering, a Gaussianregression training method is used to obtain the corresponding Gaussianprocess models and the predicted value y_(*i) of the prediction pointx_(*) under the M local Gaussian process models, and the predictionpoint generated by the clustering is subjected to weighted fusion by theposterior probability r_(*i) generated by the i^(th) Gaussian componentso as to calculate the final predicted value y_(*):

$\begin{matrix}{y_{*} = {\sum\limits_{i = 1}^{M}{r_{*i}{y_{*i}.}}}} & (18)\end{matrix}$

Further, a performance index in the following form is constructed tosolve the control input of the system:

$\begin{matrix}\begin{matrix}\min\limits_{u} & {J = {{\frac{1}{2}{\sum\limits_{i = 1}^{N}{{\xi_{k + i} - \xi_{k + i}^{d}}}_{Q_{t}}^{2}}} + {u_{k + i - 1}}_{Q_{r}}^{2}}} \\{s.t.} & {X_{k + 1} = {X_{k} + {\Delta{t \cdot {f\left( {X_{k},u_{k}} \right)}}}}} \\ & {\xi_{k} = {g\left( Q_{k} \right)}} \\ & {X_{\min} < X_{k} < X_{\max}} \\ & {\xi_{\min} < \xi_{k} < \xi_{\max}} \\ & {u_{\min} < u_{k} < u_{\max}}\end{matrix} & (19)\end{matrix}$

where Q_(r) and Q_(t) represent weighting matrices of a control errorand an input penalty term, respectively; and

-   -   by selecting appropriate controller parameters which include Δt,        N, Q_(r) and Q_(t), the above optimization control problem is        solved, the control variable U=[u_(k), u_(k+1), . . . ,        u_(k+N−1)]^(T) at the next N moments is obtained, and only the        control variable u_(k) at the moment k is used as the system        input.

Further, a generalized thruster installation matrix is given firstly,that is, there are l jet thrusters installed on the satellite platform,where the installation position of each of the thrusters is divided into[x_(Ti), y_(Ti), z_(Ti)]^(T), and the included angles between thrustvectoring and a +x direction and between the thrust vectoring and a x-yplane are respectively β_(i) and γ_(i) so that a configuration matrix ofa thruster thrust and a torque under the satellite platform system isobtained:

$\begin{matrix}{{A\left( {:,i} \right)} = \begin{bmatrix}{\cos\gamma_{i}\cos\beta_{i}} \\{\cos\gamma_{i}\sin\beta_{i}} \\{\sin\gamma_{i}} \\{{\cos\gamma_{i}\sin{\beta_{i} \cdot z_{Ti}}} - {\sin{\gamma_{i} \cdot y_{Ti}}}} \\{{\sin{\gamma_{i} \cdot x_{Ti}}} - {\cos\gamma_{i}\cos{\beta_{i} \cdot z_{Ti}}}} \\{{\cos\gamma_{i}\cos{\beta_{i} \cdot y_{Ti}}} - {\cos\gamma_{i}\sin{\beta_{i} \cdot x_{Ti}}}}\end{bmatrix}} & (20)\end{matrix}$

thus, a relationship between the platform control force and the controltorque as well as thrust of the jet thrusters is obtained:

F₀=C_(IB) AT   (21)

where C_(IB) is a conversion matrix from the platform system to aninertial system, and T=[T₁, T₂, . . . , T_(i)]^(T) is the thrust of thejet thrusters;

-   -   for the obtained platform control variable F₀, a thrust        distribution algorithm is designed by constructing the following        optimization problems, and the control variable is allocated to        each of the thrusters:

$\begin{matrix}\begin{matrix}\min\limits_{T} & {J = {\frac{1}{2}\left\lbrack {{{{C_{IB}AT} - \tau}}_{Qt}^{2} + {T}_{Rt}^{2}} \right\rbrack}} \\{s.t.} & {T_{lb} \leq T \leq T_{ub}}\end{matrix} & (22)\end{matrix}$

where T_(lb) and T_(ub) represent thrust saturation constraints for thethrusters;

-   -   the continuous thrust of each of the thrusters is obtained by        thrust distribution; in order to meet the on-off control mode of        the thrusters, the continuous thrust is converted into the        start-up time in an on-off mode through a PWM method:

$\begin{matrix}{t_{on}^{i} = {\frac{T_{i}^{\star}}{T_{i}^{o}}\Delta t}} & (23)\end{matrix}$

where T_(i)* represents the continuous thrust obtained by thedistribution algorithm, T_(i) ^(o) represents the start-up thrust of thethrusters, and Δt represents a PWM frequency; and

-   -   finally, minimum start-up time Δt_(c) of the thrusters is taken        into account, and the final thruster driving instruction is        obtained:

$\begin{matrix}{T_{on}^{i} = \left\{ {\begin{matrix}{1,} & {{{if}0} \leq t \leq {t_{on}^{i}\ {and}\ t_{on}^{i}} \geq {\Delta t_{c}}} \\{0,} & {{otherwi}se}\end{matrix}.} \right.} & (24)\end{matrix}$

The present disclosure provides electronic equipment, including a memoryand a processor, where the memory stores a computer program, and theprocessor implements the steps of the multi-arm spacecraft modelpredictive control method based on the mixture of Gaussian processeswhen executing the computer program.

The present disclosure provides a computer-readable storage medium forstoring computer instructions, where the steps of the multi-armspacecraft model predictive control method based on the mixture ofGaussian processes are implemented when the computer instructions areexecuted by the processor.

The present disclosure has the following beneficial effects.

-   -   (1) The present disclosure realizes the tracking control        capability in an on-orbit service spacecraft task space. Under a        framework of model predictive control, by considering various        interference factors and actuator saturation constraints in the        actual operating environment, the tracking of the desired        trajectory by the end-effector poses of the multi-arm spacecraft        manipulators and the platform pose is implemented, and the        design method is intuitive and simple, so that the present        disclosure is suitable for the actual engineering scene.    -   (2) For the disturbance term in the prediction model, the        Gaussian processes are designed for estimation compensation to        reduce the requirement of offline training for the amount of        data while improving the predictive control performance of the        model in the presence of various disturbances. In addition, the        Gaussian mixture model is used for clustering before training so        as to further accelerate training speed, so that the present        disclosure is more practical.    -   (3) Considering the maneuverability requirements of the on-orbit        service spacecraft during the execution of tasks, the jet        thrusters are used to control the platform pose; and since the        saturation and dead zone characteristics of the thrusters are        taken into account, the continuous control instruction generated        by the MPC is assigned as the start-up time of each of the        thrusters, thus achieving the flexibility and mobility of the        spacecraft platform.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram of a multi-arm spacecraft; and in thefigure, the multi-arm spacecraft has n manipulators, each of which has mdegrees of freedom, and has a platform equipped with l jet thrusters,where ΣB is a platform system, ΣI is an inertial system, and Σen is aend-effector coordinate system of the n^(th) manipulator.

FIG. 2 is a control flow chart of a multi-arm spacecraft task spacecontrol method based on the mixture of Gaussian processes and modelpredictive control.

DETAILED DESCRIPTION OF THE PRESENT DISCLOSURE

The technical solutions in the examples of the present disclosure willbe clearly and completely described below in conjunction with theaccompanying drawings in the examples of the present disclosure.Apparently, the described examples are only a part rather than all ofthe examples of the present disclosure. Based on the examples of thepresent disclosure, all other examples obtained by those of ordinaryskill in the art without making creative efforts shall fall within theprotection scope of the present disclosure.

Aiming at the deficiencies of the existing multi-arm on-orbit servicespacecraft control methods, the present disclosure proposes a multi-armspacecraft task space control method based on the mixture of Gaussianprocesses and model predictive control (MGP-MPC). The model predictivecontrol has excellent performance in dealing with complex nonlinearsystems such as multi-arm spacecrafts with various constraints, and iswidely applied to ground robots, unmanned aerial vehicles, autonomousdriving and other practical scenarios. Therefore, a task spacecontroller is designed based on the model predictive control. Besides,in order to enhance the anti-interference capability of the presentdisclosure, an interference model is established and compensation iscarried out in the model predictive control by utilizing thecharacteristics of small training data volume and high training speed inthe mixture of Gaussian processes. Finally, a thrust distribution methodis designed to complete platform control. The method provided by thepresent disclosure is convenient and intuitive in design and hasrelatively high practicability.

Referring to FIG. 1 and FIG. 2 , the present disclosure proposes amulti-arm spacecraft model predictive control method based on themixture of Gaussian processes, which specifically includes:

-   -   (1) establishing dynamic and kinematic models of a multi-arm        spacecraft as a prediction model in Model Predictive Control        (MPC),    -   where for a multi-arm spacecraft with n manipulators, each of        which has m degrees of freedom, and with a platform equipped        with l jet thrusters, the dynamic model of the multi-arm        spacecraft is first described in a joint space:    -   the dynamic model of the multi-arm spacecraft is expressed as:

M{umlaut over (Q)}+C+d=u   (1)

where M represents an inertial parameter of the multi-arm spacecraft; Crepresents a nonlinear term; d is the disturbance term; Q=[r₀ ^(T),q^(T)]^(T) represents information about the platform pose and jointangle, r₀ is the 6-degree-of-freedom platform pose, and q=[q₁, . . .q_(n)]^(T) is the joint angle; and u=[F₀ ^(T), F_(q) ^(T)]^(T)represents a control force and a control torque, F₀ represents aplatform control force/torque, and F_(q) represents a joint torque.

Based on a mapping relationship between the end-effectors of themanipulators and the states of joints and a platform, the descriptionform of the end-effectors poses of the multi-arm spacecraft is obtainedby using an improved DH parameter method:

$\begin{matrix}{\begin{bmatrix}r_{0} \\r_{e1} \\ \vdots \\r_{en}\end{bmatrix} = {\begin{bmatrix}r_{0} \\{{\,_{}^{1}T_{0}^{1}}{\,_{}^{1}T_{1}^{2}}\ldots{\,_{}^{1}T_{m - 1}^{m}}} \\ \vdots \\{{\,_{}^{n}T_{0}^{1}}{\,_{}^{n}T_{1}^{2}}\ldots{\,_{}^{n}T_{m - 1}^{m}}}\end{bmatrix} = {g(Q)}}} & (2)\end{matrix}$

where a homogeneous transformation matrix ^(n)T_(m-1) ^(m) is definedbased on DH parameters θ_(m),d_(m),a_(m),α_(m):

$\begin{matrix}{{\,_{}^{n}T_{m - 1}^{m}} = \begin{bmatrix}{\cos\theta_{m}} & {{- \sin}\theta_{m}} & 0 & a_{m} \\{\sin\theta_{m}\cos\alpha_{m}} & {\cos\theta_{m}\cos\alpha_{m}} & {{- \sin}\alpha_{m}} & {{- d_{i}}\sin\alpha_{m}} \\{\sin\theta_{m}\sin\alpha_{m}} & {\cos\theta_{m}\sin\alpha_{m}} & {\cos\alpha_{m}} & {d_{i}\cos\alpha_{m}} \\0 & 0 & 0 & 1\end{bmatrix}} & (3)\end{matrix}$

finally, a state space expression is built:

$\begin{matrix}\left\{ \begin{matrix}{\overset{˙}{X} = {{f\left( {X,u} \right)} = \begin{bmatrix}\overset{.}{Q} \\{M^{- 1}\left( {\tau - C - \hat{d}} \right)}\end{bmatrix}}} \\{\xi = {g(Q)}}\end{matrix} \right. & (4)\end{matrix}$

where X=[Q^(T), {dot over (Q)}^(T)]^(T), ξ=[r₀ ^(T), r_(el) ^(T), . . ., r_(en) ^(T)]^(T) represents the platform pose and the end-effectorsposes of the manipulators, and {circumflex over (d)} represents theestimated value of the disturbance term; and

-   -   sampling time Δt is set, the continuous state space expression        is converted into a discrete time model, and a system state at a        future moment k+N is predicted based on a value measured at a        current moment k:

$\begin{matrix}\left\{ {\begin{matrix}{X_{k + 1} = {X_{k} + {\Delta{t \cdot {f\left( {X_{k},u_{k}} \right)}}}}} \\{\xi_{k} = {g\left( Q_{k} \right)}}\end{matrix}.} \right. & (5)\end{matrix}$

-   -   (2) For an estimated value {circumflex over (d)} of a        disturbance term in the prediction model, inputting initial        excitation to establish a data set, and employing the mixture of        Gaussian processes for training to obtain residual dynamics of        the multi-arm spacecraft as the estimated value of the        disturbance term,    -   where the recorded p groups of data sets are first used for        training of the mixture of Gaussian processes, where an input x        contains the multi-arm spacecraft pose Q, speed {dot over (Q)},        and the increment Δu of a control variable, and an output y is a        difference value between an actual angular velocity and an        angular acceleration calculated by a nominal model (that is,        without the disturbance term), that is, the residual dynamics,        which is used to capture the impact {circumflex over (d)} of the        total disturbance on the system;    -   the Gaussian processes assume in advance that samples obey a        Gaussian distribution:

y(x)˜N (μ(x), k (x,x′))   (6)

where μ(x) represents a mean function; k (x, x′) represents a covariancefunction, which has more forms and is generally a Gaussian kernelfunction, and an expression thereof is as follows:

$\begin{matrix}{{k\left( {x,x^{\prime}} \right)} = {{\sigma_{f}^{2}{\exp\left\lbrack \frac{- \left( {x - x^{\prime}} \right)^{2}}{2l^{2}} \right\rbrack}} + {\sigma_{n}^{2}{\delta\left( {x - x^{\prime}} \right)}}}} & (7)\end{matrix}$

where σ_(n) is equal to 0.3, σ_(f) and l are hyperparameters to besolved, and

${\delta\left( {x - x^{\prime}} \right)} = \left\{ {\begin{matrix}{1,\ {x = x^{\prime}}} \\{0,\ {x \neq x^{\prime}}}\end{matrix};} \right.$

andtherefore, the Gaussian processes are actually processes of solving thehyperparameters in the covariance function through sample data. Thepresent disclosure adopts an optimal maximum likelihood estimationmethod to solve the covariance function:

$\begin{matrix}{{\max\limits_{\sigma_{f},l}{\log\left( {y❘x} \right)}} = {{{- \frac{1}{2}}y^{T}K^{- 1}y} - {\frac{1}{2}\log{❘K❘}} - {\frac{n}{2}{\log\left( {2\pi} \right)}}}} & (8)\end{matrix}$

after the hyperparameters of the Gaussian processes are obtained viatraining, the corresponding output value y, can be solved by using theinput information x, of a prediction point; a joint probability densityfunction of the samples and a test set is firstly determined:

$\begin{matrix}{\begin{bmatrix}{y(x)} \\y_{*}\end{bmatrix} \sim {N\left( {0,\begin{bmatrix}K & K_{*}^{T} \\K_{*} & K_{**}\end{bmatrix}} \right)}} & (9)\end{matrix}$ where $\begin{matrix}{K_{**} = {k\left( {x_{*},x_{*}} \right)}} & (10)\end{matrix}$ $\begin{matrix}{K_{*} = \begin{bmatrix}{k\left( {x_{*},x_{1}} \right)} & {k\left( {x_{*},x_{2}} \right)} & \ldots & {k\left( {x_{*},x_{p}} \right)}\end{bmatrix}} & (11)\end{matrix}$ $\begin{matrix}{K = \begin{bmatrix}{k\left( {x_{1},x_{1}} \right)} & {k\left( {x_{1},x_{2}} \right)} & \ldots & {k\left( {x_{1},x_{p}} \right)} \\{k\left( {x_{2},x_{1}} \right)} & {k\left( {x_{2},x_{2}} \right)} & \ldots & {k\left( {x_{2},x_{p}} \right)} \\ \vdots & \vdots & \ddots & \vdots \\{k\left( {x_{p},x_{1}} \right)} & {k\left( {x_{p},x_{2}} \right)} & \ldots & {k\left( {x_{p},x_{p}} \right)}\end{bmatrix}} & (12)\end{matrix}$

-   -   a posterior probability expression of the prediction point is        obtained based on Bayesian reasoning:

y _(*) |y(x)˜N(K _(*) K ⁻¹ y, K _(**) −K _(*) K ⁻¹ K _(*) ^(T))   (13)

a predicted value of the Gaussian processes is accordingly obtained:

y _(*) =K _(*) K ⁻¹ y   (14).

On the basis of the above-mentioned Gaussian processes of the multi-armspacecraft, in order to speed up training and predicting, the presentdisclosure utilizes a local Gaussian process method in the Gaussianprocesses of the multi-arm spacecraft. The local Gaussian process methodincludes: clustering based on a Gaussian mixture model method,performing regression on each of Gaussian components, and finally fusinglocal predicted values during prediction;

-   -   assuming that the input x obeys the mixture distribution of the        M Gaussian components, a probability density distribution        function of the input x is as follows:

$\begin{matrix}{{p\left( {x❘\Theta} \right)} = {\sum\limits_{i = 1}^{M}{\pi_{i}{p\left( {{x❘\mu_{i}},\sum_{i}} \right)}}}} & (15)\end{matrix}$

where π_(i) is a mixing coefficient, which represents the probabilitythat the sample comes from the i^(th) Gaussian component; p(x|μ_(i),Σ_(i)) represents the probability that the sample is generated by thei^(th) Gaussian component; and {μ_(i), Σ_(i)} is the hyperparameter ofthe mixture of Gaussian processes.

EM (Expectation-Maximization) is an effective algorithm for learningmixture models, which is theoretically equivalent to a maximumlikelihood estimation method. After the training samples are obtained,the model parameters are first initialized, and then the modelparameters are updated by means of the continuous iteration of steps Eand M of expectation-maximizationuntil the model converges,

-   -   where step E uses the parameters {π_(i), μ_(i), Σ_(i)} estimated        in step M to calculate a posterior probability r_(i) of the        sample coming from the i^(th) Gaussian component:

$\begin{matrix}{r_{i} = \frac{\pi_{i}{p\left( {{x❘\mu_{i}},\sum_{i}} \right)}}{\sum\limits_{j = 1}^{M}{\pi_{j}{p\left( {{x❘\mu_{j}},\sum_{j}} \right)}}}} & (16)\end{matrix}$

in step M, the posterior probability r_(i) of the training samples iscalculated through step E, and the model parameters are updated bymaximizing a logarithmic likelihood function:

$\begin{matrix}\left\{ \begin{matrix}{\mu_{i} = {\sum\limits_{k = 1}^{n}{r_{k}{x_{k}/{\sum\limits_{k = 1}^{n}r_{k}}}}}} \\{\sum_{i}{= {\sum\limits_{k = 1}^{n}{{r_{k}\left( {x_{k} - \mu_{i}} \right)}{\left( {x_{k} - \mu_{i}} \right)^{T}/{\sum\limits_{k = 1}^{n}r_{k}}}}}}} \\{\pi_{i} = {\frac{1}{n}{\sum\limits_{k = 1}^{n}r_{k}}}}\end{matrix} \right. & (17)\end{matrix}$

for each group of Gaussian components obtained by clustering, a Gaussianregression training method is used to obtain the corresponding Gaussianprocess models and the predicted value y_(*i) of the prediction pointx_(*) under the M local Gaussian process models, and the predictionpoint generated by the clustering is subjected to weighted fusion by theposterior probability r_(*i) generated by the i^(th) Gaussian componentso as to calculate the final predicted value y_(*):

$\begin{matrix}{y_{*} = {\sum\limits_{i = 1}^{M}{r_{*i}{y_{*i}.}}}} & (18)\end{matrix}$

-   -   (3) Constructing, according to actual drive input saturation        constraints, a model predictive controller to implement the        tracking of a desired trajectory by end-effectors poses of        multi-arm spacecraft manipulators and a platform pose, where the        saturation includes saturation of joint motors and thrust        saturation of platform thrusters,    -   in order to reduce the energy consumption as much as possible        while realizing the task space trajectory tracking of the        on-orbit service spacecraft, a performance index in the        following form is constructed to solve the control input of the        system:

$\begin{matrix}\begin{matrix}\min\limits_{u} & {J = {{\frac{1}{2}{\sum\limits_{i = 1}^{N}{{\xi_{k + i} - \xi_{k + i}^{d}}}_{Q_{t}}^{2}}} + {u_{k + i - 1}}_{Q_{r}}^{2}}}\end{matrix} & (19)\end{matrix}$ $\begin{matrix}{s.t.} & {X_{k + 1} = {X_{k} + {{\Delta t} \cdot {f\left( {X_{k},\ u_{k}} \right)}}}}\end{matrix}$ ξ_(k) = g(Q_(k)) X_(min) < X_(k) < X_(max)ξ_(min) < ξ_(k) < ξ_(max) u_(min) < u_(k) < u_(max)

where Q_(r) and Q_(t) represent weighting matrices of a control errorand an input penalty term, respectively; and

-   -   by selecting appropriate controller parameters which include Δt,        N, Q_(r) and Q_(t), the above optimization control problem is        solved, the control variable U=[u_(k), u_(k+1), . . . ,        u_(k+N−1)]^(T) at the next N moments is obtained, and only the        control variable u_(k) at the moment k is used as the system        input.    -   (4) Setting the platform thrusters to be in an on-off control        mode, and assigning a platform continuous control instruction        generated by the M PC as start-up time of each of the thrusters        to obtain a final thruster driving instruction.

A generalized thruster installation matrix is given first, that is,there are l jet thrusters installed on the satellite platform, where theinstallation position of each of the thrusters is divided into [x_(Ti),y_(Ti), z_(Ti)]^(T), and the included angles between thrust vectoringand a +x direction and between the thrust vectoring and a x-y plane arerespectively β_(i) and γ_(i) so that a configuration matrix of athruster thrust and a torque under the satellite platform system isobtained:

$\begin{matrix}{{A\left( {:,i} \right)} = \begin{bmatrix}{\cos\gamma_{i}\cos\beta_{i}} \\{\cos\gamma_{i}\sin\beta_{i}} \\{\sin\gamma_{i}} \\{{\cos\gamma_{i}\sin{\beta_{i} \cdot z_{Ti}}} - {\sin{\gamma_{i} \cdot y_{Ti}}}} \\{{\sin{\gamma_{i} \cdot x_{Ti}}} - {\cos\gamma_{i}\cos{\beta_{i} \cdot z_{Ti}}}} \\{{\cos\gamma_{i}\cos{\beta_{i} \cdot y_{Ti}}} - {\cos\gamma_{i}\sin{\beta_{i} \cdot x_{Ti}}}}\end{bmatrix}} & (20)\end{matrix}$

thus, a relationship between the platform control force and the controltorque as well as thrust of the jet thrusters is obtained:

F₀=C_(IB)AT   (21)

where C_(IB) is a conversion matrix from the platform system to aninertial system, and T=[T₁, T₂, . . . , T_(i)]^(T) is the thrust of thejet thrusters;

-   -   for the obtained platform control variable F₀, a thrust        distribution algorithm is designed by constructing the following        optimization problems, and the control variable is allocated to        each of the thrusters:

$\begin{matrix}{\begin{matrix}\min\limits_{T} & {J = \frac{1}{2}}\end{matrix}\left\lbrack {{{{C_{IB}{AT}} - \tau}}_{Qt}^{2} + {T}_{Rt}^{2}} \right\rbrack} & (22)\end{matrix}$ $\begin{matrix}{s.t.} & {T_{lb} \leq T \leq T_{ub}}\end{matrix}$

where T_(lb) and T_(ub) represent thrust saturation constraints for thethrusters;

-   -   the continuous thrust of each of the thrusters is obtained by        thrust distribution; in order to meet the on-off control mode of        the thrusters, the continuous thrust is converted into the        start-up time in an on-off mode through a PWM method:

$\begin{matrix}{t_{on}^{i} = {\frac{T_{i}^{\star}}{T_{i}^{o}}\Delta t}} & (23)\end{matrix}$

where T_(i)* represents the continuous thrust obtained by thedistribution algorithm, T_(i) ^(o) represents the start-up thrust of thethrusters, and Δt represents a PWM frequency; and

-   -   finally, minimum start-up time Δt_(c) of the thrusters is taken        into account, and the final thruster driving instruction is        obtained:

$\begin{matrix}{T_{on}^{i} = \left\{ {\begin{matrix}{1,} & {{{if}0} \leq t \leq {t_{on}^{i}{and}t_{on}^{i}} \geq {\Delta t_{c}}} \\{0,} & {otherwise}\end{matrix}.} \right.} & (24)\end{matrix}$

The present disclosure provides electronic equipment, including a memoryand a processor, where the memory stores a computer program, and theprocessor implements the steps of the multi-arm spacecraft modelpredictive control method based on the mixture of Gaussian processeswhen executing the computer program.

The present disclosure provides a computer-readable storage medium forstoring computer instructions, where the steps of the multi-armspacecraft model predictive control method based on the mixture ofGaussian processes are implemented when the computer instructions areexecuted by the processor.

The multi-arm spacecraft model predictive control method based on themixture of Gaussian processes, the equipment, and the medium which areproposed by the present disclosure have been described above in detail.Specific examples are used herein to illustrate the principle andimplementation of the present disclosure. The description of the aboveexamples is only used to help understand the method and its core idea ofthe present disclosure. Furthermore, for those skilled in the art,according to the idea of the present disclosure, there may be changes inthe specific examples and the scope of application. To sum up, thecontent of the Description should not be construed as limiting thepresent disclosure.

1. A multi-arm spacecraft model predictive control method based on amixture of Gaussian processes, comprising: establishing dynamic andkinematic models of a multi-arm spacecraft as a prediction model in MPC;for an estimated value {circumflex over (d)}of a disturbance term in theprediction model, inputting initial excitation to establish a data set,and employing the mixture of Gaussian processes for training to obtainresidual dynamics of the multi-arm spacecraft as the estimated value ofthe disturbance term; constructing, according to actual drive inputsaturation constraints, a model predictive controller to implementtracking of a desired trajectory by end-effectors poses of spacecraftmanipulators and a platform pose, wherein the saturation comprisessaturation of joint motors and thrust saturation of platform thrusters;and setting the platform thrusters to be in an on-off control mode, andassigning a platform continuous control instruction generated by the MPCas start-up time of each of the thrusters to obtain a final thrusterdriving instruction.
 2. The method according to claim 1, wherein for amulti-arm spacecraft with n manipulators, each of which has m degrees offreedom, and with a platform equipped with l jet thrusters, the dynamicmodel of the multi-arm spacecraft is first described in a joint space:the dynamic model of the multi-arm spacecraft is expressed as:M{umlaut over (Q)}+C+d=u   (1), where M represents an inertial parameterof the multi-arm spacecraft; C represents a nonlinear term; d is thedisturbance term; Q=[r₀ ^(T), q^(T)]^(T) represents information aboutthe platform pose and joint angle, r₀ is a 6-degree-of-freedom platformpose, and q=[q₁, . . . , q_(n)]^(T) is the joint angle; and u=[F₀ ^(T),F_(q) ^(T)]^(T) represents a control force and a control torque, F₀represents a platform control force/torque, and F_(q) represents a jointtorque.
 3. The method according to claim 2, wherein based on a mappingrelationship between the end-effectors of the manipulators and states ofjoints and a platform, a description form of the end-effectors poses ofthe multi-arm spacecraft manipulators is obtained by using a DHparameter method: [ r 0 r e ⁢ 1 ⋮ r en ] = [ r 0 1 T 0 1 1 T 1 2 ⁢ … 1 Tm - 1 m ⋮ n T 0 1 n T 1 2 ⁢ … n T m - 1 m ] = g ⁡ ( Q ) , ( 2 ) where ahomogeneous transformation matrix ^(n)T_(m-1) ^(m) is defined based onDH parameters θ_(m),d_(m),a_(m),α_(m): n T m - 1 m = [ cos ⁢ θ m - s ⁢ in ⁢θ m 0 a m sin ⁢ θ m ⁢ cos ⁢ α m cos ⁢ θ m ⁢ cos ⁢ α m - s ⁢ in ⁢ α m - d i ⁢ sin ⁢α m sin ⁢ θ m ⁢ sin ⁢ α m cos ⁢ θ m ⁢ sin ⁢ α m cos ⁢ α m d i ⁢ cos ⁢ α m 0 0 0 1] , ( 3 ) finally, a state space expression is built: $\begin{matrix}\left\{ {\begin{matrix}{\overset{.}{X} = {{f\left( {X,u} \right)} = \begin{bmatrix}\overset{.}{Q} \\{M^{- 1}\left( {\tau - C - \overset{\bigwedge}{\left. d \right)}} \right.}\end{bmatrix}}} \\{\xi = {g(Q)}}\end{matrix},} \right. & (4)\end{matrix}$ where X=[Q^(T), {dot over (Q)}^(T)]^(T), ξ=[r₀ ^(T),r_(el) ^(T), . . . , r_(en) ^(T)]^(T) represents the platform pose andthe end-effectors poses of the manipulators, and {circumflex over (d)}represents the estimated value of the disturbance term; and samplingtime Δt is set, continuous state space expression is converted into adiscrete time model, and a system state at a future moment k+N ispredicted based on a value measured at a current moment k:$\begin{matrix}\left\{ {\begin{matrix}{X_{k + 1} = {X_{k} + {\Delta{t \cdot {f\left( {X_{k},u_{k}} \right)}}}}} \\{\xi = {g\left( Q_{k} \right)}}\end{matrix}.} \right. & (5)\end{matrix}$
 4. The method according to claim 3, wherein recorded pgroups of data sets are used for training of the mixture of Gaussianprocesses, wherein an input x contains the multi-arm spacecraft pose Q,speed {dot over (Q)}, and an increment Δu of a control variable, and anoutput y is a difference value between an actual angular velocity and anangular acceleration calculated by a nominal model, that is, theresidual dynamics, which is used to capture the impact {circumflex over(d)} of a total disturbance on the system; the Gaussian processes assumein advance that samples obey a Gaussian distribution:y(x)˜N (μ(x), k (x,x′))   (6), where μ(x) represents a mean function; k(x , x′) represents a covariance function, and an expression thereof isas follows: $\begin{matrix}{{{k\left( {x,x^{\prime}} \right)} = {{\sigma_{f}^{2}{\exp\left\lbrack \frac{- \left( {x - x^{\prime}} \right)^{2}}{2l^{2}} \right\rbrack}} + {\sigma_{n}^{2}{\delta\left( {x - x^{\prime}} \right)}}}},} & (7)\end{matrix}$ where σ_(n) is equal to 0.3, σ_(f) and l arehyperparameters to be solved, and${\delta\left( {x - x^{\prime}} \right)} = \left\{ {\begin{matrix}{1,\ {x = x^{\prime}}} \\{0,\ {x \neq x^{\prime}}}\end{matrix};} \right.$ the covariance function is solved by using anoptimal maximum likelihood estimation method: $\begin{matrix}{{{\max\limits_{\sigma_{f},l}{\log\left( y \middle| x \right)}} = {{{- \frac{1}{2}}y^{T}K^{1}y} - {\frac{1}{2}\log{❘K❘}} - {\frac{n}{2}{\log\left( {2\pi} \right)}}}},} & (8)\end{matrix}$ after the hyperparameters of the Gaussian processes areobtained via training, the corresponding output value y_(*) can besolved by using the input information x_(*) of a prediction point; ajoint probability density function of the samples and a test set isfirstly determined: $\begin{matrix}{{\begin{bmatrix}{y(x)} \\y_{*}\end{bmatrix} \sim {N\left( {\theta,\begin{bmatrix}K & K_{\star}^{T} \\K_{*} & K_{**}\end{bmatrix}} \right)}},{where}} & (9) \\{K_{\star *} = {k\left( {x_{\star},x_{\star}} \right)}} & (10) \\{K_{*} = \left\lbrack {{k\left( {x_{*},x_{1}} \right)}\ {k\left( {x_{*},x_{2}} \right)}\ {k\left( {x_{*},x_{p}} \right)}} \right\rbrack} & (11)\end{matrix}$ a posterior probability expression of the prediction pointis obtained based on Bayesian reasoning:y _(*) |y(x)˜N(K _(*) K ⁻¹ y, K _(**) −K _(*) K ⁻¹ K _(*) ^(T))   (13),a predicted value of the Gaussian processes is accordingly obtained:y _(*) =K _(*) K ⁻¹ y   (14).
 5. The method according to claim 4,wherein a local Gaussian process method is utilized in the Gaussianprocesses of the multi-arm spacecraft, and the local Gaussian processmethod comprises: clustering based on a Gaussian mixture model method,performing regression on each of Gaussian components, and finally fusinglocal predicted values during prediction; assuming that the input xobeys mixture distribution of M Gaussian components, a probabilitydensity distribution function of the input x is as follows:$\begin{matrix}{{p\left( {x❘\Theta} \right)} = {\overset{M}{\sum\limits_{i = 1}}{\pi_{i}{p\left( {{x❘\mu_{i}},\Sigma_{i}} \right)}}}} & (15)\end{matrix}$ where π_(i) is a mixing coefficient, which represents theprobability that the sample comes from the i^(th) Gaussian component;p(x|μ_(i), Σ_(i)) represents the probability that the sample isgenerated by the i^(th) Gaussian component; and {μ_(i), Σ_(i)} is thehyperparameter of the mixture of Gaussian processes.
 6. The methodaccording to claim 5, wherein after the training samples are obtained,model parameters are first initialized, and then the model parametersare updated by means of continuous iteration of steps E and M ofexpectation-maximization until the model converges; where step E usesthe parameters {π_(i), μ_(i), Σ_(i)} estimated in step M to calculate aposterior probability r_(i) of the sample coming from the i^(th)Gaussian component: $\begin{matrix}{{r_{i} = \frac{\pi_{i}{p\left( {{x❘\mu_{i}},\Sigma_{i}} \right)}}{\overset{M}{\sum\limits_{j = 1}}{\pi_{i}{p\left( {{x❘\mu_{j}},\Sigma_{j}} \right)}}}},} & (16)\end{matrix}$ in step M, the posterior probability r_(i) of the trainingsamples is calculated through step E, and the model parameters areupdated by maximizing a logarithmic likelihood function: $\begin{matrix}\left\{ {\begin{matrix}{\mu_{i} = {\overset{n}{\sum\limits_{k = 1}}{r_{k}{x_{k}/{\overset{n}{\sum\limits_{k = 1}}r_{k}}}}}} \\{\Sigma_{i} = {\overset{n}{\sum\limits_{k = 1}}{{r_{k}\left( {x_{k} - \mu_{i}} \right)}{\left( {x_{k} - \mu_{i}} \right)^{T}/{\overset{n}{\sum\limits_{k = 1}}r_{k}}}}}} \\{\pi_{i} = {\frac{1}{n}{\overset{n}{\sum\limits_{k = 1}}r_{k}}}}\end{matrix},} \right. & (17)\end{matrix}$ for each group of Gaussian components obtained byclustering, a Gaussian regression training method is used to obtain thecorresponding Gaussian process models and the predicted value y_(*i) ofthe prediction point x_(*) under the M local Gaussian process models,and the prediction point generated by the clustering is subjected toweighted fusion by the posterior probability r_(*i) generated by thei^(th) Gaussian component so as to calculate the final predicted valuey_(*): $\begin{matrix}{y_{*} = {\overset{M}{\sum\limits_{i = 1}}{r_{*i}{y_{*i}.}}}} & (18)\end{matrix}$
 7. The method according to claim 6, wherein a performanceindex in the following form is constructed to solve the control input ofthe system: $\begin{matrix}{{\min\limits_{u}J} = {{\frac{1}{2}{\sum\limits_{i = 1}^{N}{\|\xi_{k + i}}}} - {\xi_{k + i}^{d}\|_{Q_{t}}^{2}} + {u_{k + i - 1}}_{Q_{r}}^{2}}} & (19)\end{matrix}$ s.t.X_(k + 1) = X_(k) + Δt ⋅ f(X_(k), u_(k))ξ_(k) = g(Q_(k)) X_(min) < X_(k) < X_(max) ξ_(min) < ξ_(k) < ξ_(max)u_(min) < u_(k) < u_(max), where Q_(r) and Q_(t) represent weightingmatrices of a control error and an input penalty term, respectively; andby selecting appropriate controller parameters which comprise Δt, N,Q_(r) and Q_(t), the above optimization control problem is solved, thecontrol variable U=[u _(k), u_(k+1), . . . , u_(k+N−1)]^(T) at the nextN moments is obtained, and only the control variable u_(k) at the momentk is used as the system input.
 8. The method according to claim 7,wherein a generalized thruster installation matrix is given first, thatis, there are l jet thrusters installed on the satellite platform,wherein the installation position of each of the thrusters is dividedinto [x _(Ti), y_(Ti), z_(Ti)]^(T), and the included angles betweenthrust vectoring and a +x direction and between the thrust vectoring anda x-y plane are respectively β_(i) and γ_(i) so that a configurationmatrix of a thruster thrust and a torque under the satellite platformsystem is obtained: $\begin{matrix}{{{A\left( {:{,i}} \right)} = \begin{bmatrix}{\cos\gamma_{i}\cos\beta_{i}} \\{\cos\gamma_{i}\sin\beta_{i}} \\{\sin\gamma_{i}} \\{{\cos\gamma_{i}\sin{\beta_{i} \cdot z_{Ti}}} - {\sin{\gamma_{i} \cdot y_{Ti}}}} \\{{\sin{\gamma_{i} \cdot z_{Ti}}} - {\cos\gamma_{i}\cos{\beta_{i} \cdot z_{Ti}}}} \\{{\cos\gamma_{i}\cos{\beta_{i} \cdot y_{Ti}}} - {\cos\gamma_{i}\sin{\beta_{i} \cdot x_{Ti}}}}\end{bmatrix}},} & (20)\end{matrix}$ thus, a relationship between the platform control forceand the control torque as well as thrust of the jet thrusters isobtained:F₀=C_(IB)AT   (21), where C_(IB) is a conversion matrix from theplatform system to an inertial system, and T=[T₁, T₂, . . . , T_(l)]^(T)is the thrust of the jet thrusters; or the obtained platform controlvariable F₀, a thrust distribution algorithm is designed by constructingthe following optimization problems, and the control variable isallocated to each of the thrusters: $\begin{matrix}{{{\min\limits_{T}J} = {\frac{1}{2}\left\lbrack {{{{C_{IB}AT} - \tau}}_{Qt}^{2} + {T}_{Rf}^{2}} \right\rbrack}}{{{s.t.T_{lb}} \leq T \leq T_{ub}},}} & (22)\end{matrix}$ where T_(lb) and T_(ub) represent thrust saturationconstraints for the thrusters; the continuous thrust of each of thethrusters is obtained by thrust distribution; in order to meet theon-off control mode of the thrusters, the continuous thrust is convertedinto the start-up time in an on-off mode through a PWM method:$\begin{matrix}{{t_{on}^{i} = {\frac{T_{i}^{*}}{T_{i}^{o}}\Delta t}},} & (23)\end{matrix}$ where T_(i)* represents the continuous thrust obtained bythe distribution algorithm, T_(i) ^(o) represents the start-up thrust ofthe thrusters, and Δt represents a PWM frequency; and finally, minimumstart-up time Δt_(c) of the thrusters is taken into account, and thefinal thruster driving instruction is obtained: $\begin{matrix}{T_{on}^{i} = \left\{ {\begin{matrix}{1,\ {{{if}\ 0} \leq t \leq {t_{on}^{i}\ {and}\ t_{on}^{i}} \geq {\Delta t_{c}}}} \\{0,\ {otherwise}}\end{matrix}.} \right.} & (24)\end{matrix}$
 9. Electronic equipment, comprising a memory and aprocessor, wherein the memory stores a computer program, and theprocessor implements the steps of the method according to claim 1 whenexecuting the computer program.
 10. A computer-readable storage mediumfor storing computer instructions, wherein the steps of the methodaccording to claim 1 are implemented when the computer instructions areexecuted by the processor.